##########################################
## "Voters get what they want"          ##
##########################################
## Heinrich, Kobayashi, Long            ##
##########################################

##          
## Re-analyze Lake/Nielsen data 
## Script 02
## June 29, 2017

## Differences from Nielsen: Only USA, data ends in 2000,
## recipient random effects


## Load data and set things up
##############################
load("output/data_SH.Rdata")
load("output/data_EH.Rdata")

## Set up storages
out_lmer <- out <- vector("list", 2)
for(i in 1:2) out[[i]] <- vector("list", chains)


### Estimate all models
#######################
for(k in 1:chains)
{
  ## Total aid, SH sample
  out_tmp <- MCMCglmm(cbind(lnaidpc_l, lnaidpc) ~ 1 + lphysint + polity2 + llnaidpc + lnworldaidtotal + lln_rgdpc + lln_population +
                        + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar +  
                        + ldonorallyneighbor2 + ls3un + llnNYT_HR + llnreftotal,
                      random= ~ us(1):countrynumcode_g,
                      data=data_SH, verbose=round(niter/5), 
                      family="cengaussian", nitt=niter, burnin=burnin, thin=1)
  out[[1]][[k]] <- mcmc(cbind(out_tmp$Sol, out_tmp$VCV))
  out_lmer[[1]] <- lmer(lnaidpc ~ (1 | countrynumcode_g) + lphysint + polity2 + llnaidpc + lnworldaidtotal + lln_rgdpc + lln_population +
                          + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar + 
                          + ldonorallyneighbor2 + ls3un + llnNYT_HR + llnreftotal,
                        data=data_SH)
  
  ## TotalAid, EH sample
  out_tmp <- MCMCglmm(cbind(lnaidpc_l, lnaidpc) ~ 1 + lphysint + polity2 + llnaidpc + lnworldaidtotal + lln_rgdpc + lln_population# +
                      + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar + 
                        + ldonorallyneighbor2 + ls3un + llnNYT_HR + llnreftotal,
                      random= ~ us(1):countrynumcode_g,
                      prior=list(R=list(V=1, nu=.001), G=list(G1=list(V=1, nu=.001))),
                      data=data_EH, verbose=round(niter/5), 
                      family="cengaussian", nitt=niter, burnin=burnin, thin=1)
  out[[2]][[k]] <- mcmc(cbind(out_tmp$Sol, out_tmp$VCV))
  out_lmer[[2]] <- lmer(lnaidpc ~ (1 | countrynumcode_g) + lphysint + polity2 + llnaidpc + lnworldaidtotal + lln_rgdpc + lln_population +
                          + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar + 
                          + ldonorallyneighbor2 + ls3un + llnNYT_HR + llnreftotal,
                        data=data_EH)
}


## Convergence
##############
load("output/Rhat.Rdata")
Rhat <- rbind(Rhat, data.frame(Rhat=round(max(c(laply(.data=out, 
                                                      .fun=function(x) gelman.diag(x)$psrf[,1])), decreasing=TRUE),
                                          di=2),
                               Which="Nielsen re-analysis"))
save(Rhat, file="output/Rhat.Rdata")


## Make regression table
########################
out2 <- llply(.data=out, .fun=function(x) ldply(.data=x, .fun=function(y) y))
mods_ci <- mods_coefs <- mods_tmp <- vector("list", 2)
mods_vcv_1 <- mods_vcv_2 <- matrix("", 2, 2)

for(i in 1:2)
{
  tmp_coef <- as.matrix(out2[[i]][, 1:(ncol(out2[[i]])-2)])
  res_var <- out2[[i]][, ncol(out2[[i]])]
  re_var <- out2[[i]][, ncol(out2[[i]])-1]
  mods_tmp[[i]] <- out_lmer[[i]]
  mods_coefs[[i]] <- colMeans(tmp_coef)
  mods_ci[[i]] <- colQuantiles(tmp_coef, probs=c(.025, .975))
  mods_vcv_1[i,1] <- round(mean(res_var), di=2)
  mods_vcv_2[i,1] <- round(mean(re_var), di=2)
  mods_vcv_1[i,2] <- paste0("(", round(quantile(res_var, probs=.025,), di=2),
                            ", ", round(quantile(res_var, probs=.975,), di=2), ")")
  mods_vcv_2[i,2] <- paste0("(", round(quantile(re_var, probs=.025,), di=2),
                            ", ", round(quantile(re_var, probs=.975,), di=2), ")")
}

table_nielsen <- stargazer(type="latex", style="ajps",
                     mods_tmp, coef=mods_coefs, ci.custom=mods_ci, 
                     ci=TRUE,
                     p.auto=FALSE, report="vcs", table.placement="t",
                     dep.var.labels="Total aid (Nielsen 2013)",label="TabA1:Nielsen",
                     column.labels=c("Security hierarchy sample", "Economic hierarchy sample"),
                     font.size="tiny", model.numbers=FALSE, model.names=FALSE,
                     covariate.labels=c("Human rights violations", "Democracy", "Lagged total aid", "Lagged economic aid",  
                                        "Global aid flows", "GDP per capita", "Population", "Trade", "Alliance", "Socialist",
                                        "Cold War", "Cold War x socialist", "War", "Ally neighbor", "UN voting similarity", "News (human rights)", 
                                        "Refugees", "Constant"),
                     digits=2, notes="", df=FALSE, omit.stat=c("ll", "adj.rsq", "f", "rsq"))
table_nielsen <- table_nielsen[-5]

which_row <- which(grepl(x=table_nielsen, pattern="AIC") == TRUE)
table_nielsen <- c(table_nielsen[1:(which_row-1)],  
                   paste("Residual SE & ", paste(mods_vcv_1[,1], collapse="  & "), "\\\\", sep=""),
                   paste("& ", paste(mods_vcv_1[,2], collapse="  & "), "\\\\", sep=""),
                   "  & & \\\\ ",
                   paste("Random effect SE & ", paste(mods_vcv_2[,1], collapse="  & "), "\\\\", sep=""),
                   paste(" & ", paste(mods_vcv_2[,2], collapse="  & "), "\\\\", sep=""),
                   "\\hline \\\\[-1.8ex] ", 
                   "\\end{tabular} ",
                   "\\caption{\\textbf{Posterior summaries for the replication of Nielsen (2013) using our hierarchy subsamples}. Each column depicts 
                    the results for a different hierarchy data subset. The standalone number is the median posterior
                   estimate, the range below it the 95\\% central credible interval.}",
                   "\\end{table} ")
table_nielsen <- str_replace_all(string=table_nielsen, pattern="\\[t\\] ", replacement="\\[!ht\\] ")
tmp <- table_nielsen[5]
table_nielsen <- table_nielsen[-5]
table_nielsen[length(table_nielsen)-1] <- paste0(table_nielsen[length(table_nielsen)-1], tmp)
write(table_nielsen, file="output/tables/table_a1_nielsen.tex")
